{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8f40cb8f-9104-48ce-8ee6-fffa4671045c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "import pandas as pd\n",
    "from scipy.stats import spearmanr\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set_style('whitegrid')\n",
    "\n",
    "def p_adjust(pvalues, method='fdr_bh'):\n",
    "    res = multipletests(pvalues, method=method)\n",
    "    return np.array(res[1], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3af6be6-d931-4378-9330-a064cbc0fb22",
   "metadata": {},
   "source": [
    "# Correlate nasal 16S with vaccine response\n",
    "\n",
    "##### Michael Shaffer\n",
    "##### 7/21/22\n",
    "##### Merck ESC, Sys bio group\n",
    "\n",
    "To look for associations between the nasal microbiome and vaccine response we have calculated correlations between the abundances of individual OTUs and the continuous titer measurements from 1 year of life."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d18d0e-aa22-4f85-be4e-ed4c6f4a2516",
   "metadata": {},
   "source": [
    "## Read in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "663ff2b9-d055-41a8-80c7-70179b7c7872",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SubmissionType</th>\n",
       "      <th>SampleNumber</th>\n",
       "      <th>SampleIDValidation</th>\n",
       "      <th>DiversigenCheckInSampleName</th>\n",
       "      <th>ReplacesLowVolumeSampleID</th>\n",
       "      <th>BoxLocation</th>\n",
       "      <th>SampleType</th>\n",
       "      <th>SampleSource</th>\n",
       "      <th>SequencingType</th>\n",
       "      <th>BabyN</th>\n",
       "      <th>...</th>\n",
       "      <th>median_mmNorm_PCV</th>\n",
       "      <th>median_mmNorm_DTAPHib</th>\n",
       "      <th>protectNorm_Dip</th>\n",
       "      <th>protectNorm_TET</th>\n",
       "      <th>protectNorm_PRP (Hib)</th>\n",
       "      <th>protectNorm_PT</th>\n",
       "      <th>protectNorm_PRN</th>\n",
       "      <th>protectNorm_FHA</th>\n",
       "      <th>geommean_protectNorm</th>\n",
       "      <th>VR_group_v2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>106_V5_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A3</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>106</td>\n",
       "      <td>...</td>\n",
       "      <td>0.061955</td>\n",
       "      <td>0.052874</td>\n",
       "      <td>2.1</td>\n",
       "      <td>3.0</td>\n",
       "      <td>2.600000</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>1.3750</td>\n",
       "      <td>1.140388</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V2_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A4</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>107</td>\n",
       "      <td>...</td>\n",
       "      <td>0.958142</td>\n",
       "      <td>0.114018</td>\n",
       "      <td>4.4</td>\n",
       "      <td>5.2</td>\n",
       "      <td>10.666667</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>1.1250</td>\n",
       "      <td>0.3750</td>\n",
       "      <td>1.783418</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V3_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>107_V8_NS_A1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A5</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>107</td>\n",
       "      <td>...</td>\n",
       "      <td>0.958142</td>\n",
       "      <td>0.114018</td>\n",
       "      <td>4.4</td>\n",
       "      <td>5.2</td>\n",
       "      <td>10.666667</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>1.1250</td>\n",
       "      <td>0.3750</td>\n",
       "      <td>1.783418</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V5_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A8</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>107</td>\n",
       "      <td>...</td>\n",
       "      <td>0.958142</td>\n",
       "      <td>0.114018</td>\n",
       "      <td>4.4</td>\n",
       "      <td>5.2</td>\n",
       "      <td>10.666667</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>1.1250</td>\n",
       "      <td>0.3750</td>\n",
       "      <td>1.783418</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108_V4_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A9</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>108</td>\n",
       "      <td>...</td>\n",
       "      <td>0.003102</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.5</td>\n",
       "      <td>1.800000</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>0.3125</td>\n",
       "      <td>0.1875</td>\n",
       "      <td>0.449420</td>\n",
       "      <td>LVR</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 84 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               SubmissionType  SampleNumber  SampleIDValidation  \\\n",
       "SampleID                                                          \n",
       "106_V5_NS_A1  Primary in Tube             2                 NaN   \n",
       "107_V2_NS_A1  Primary in Tube             3                 NaN   \n",
       "107_V3_NS_A1  Primary in Tube             4                 NaN   \n",
       "107_V5_NS_A1  Primary in Tube             5                 NaN   \n",
       "108_V4_NS_A1  Primary in Tube             6                 NaN   \n",
       "\n",
       "             DiversigenCheckInSampleName ReplacesLowVolumeSampleID  \\\n",
       "SampleID                                                             \n",
       "106_V5_NS_A1                         NaN                       NaN   \n",
       "107_V2_NS_A1                         NaN                       NaN   \n",
       "107_V3_NS_A1                107_V8_NS_A1                       NaN   \n",
       "107_V5_NS_A1                         NaN                       NaN   \n",
       "108_V4_NS_A1                         NaN                       NaN   \n",
       "\n",
       "             BoxLocation  SampleType  SampleSource SequencingType  BabyN  ...  \\\n",
       "SampleID                                                                  ...   \n",
       "106_V5_NS_A1   Box 1, A3  Nasal Swab  Human Infant            16S    106  ...   \n",
       "107_V2_NS_A1   Box 1, A4  Nasal Swab  Human Infant            16S    107  ...   \n",
       "107_V3_NS_A1   Box 1, A5  Nasal Swab  Human Infant            16S    107  ...   \n",
       "107_V5_NS_A1   Box 1, A8  Nasal Swab  Human Infant            16S    107  ...   \n",
       "108_V4_NS_A1   Box 1, A9  Nasal Swab  Human Infant            16S    108  ...   \n",
       "\n",
       "              median_mmNorm_PCV median_mmNorm_DTAPHib  protectNorm_Dip  \\\n",
       "SampleID                                                                 \n",
       "106_V5_NS_A1           0.061955              0.052874              2.1   \n",
       "107_V2_NS_A1           0.958142              0.114018              4.4   \n",
       "107_V3_NS_A1           0.958142              0.114018              4.4   \n",
       "107_V5_NS_A1           0.958142              0.114018              4.4   \n",
       "108_V4_NS_A1           0.003102              0.000000              0.5   \n",
       "\n",
       "             protectNorm_TET  protectNorm_PRP (Hib) protectNorm_PT  \\\n",
       "SampleID                                                             \n",
       "106_V5_NS_A1             3.0               2.600000         0.3125   \n",
       "107_V2_NS_A1             5.2              10.666667         0.3125   \n",
       "107_V3_NS_A1             5.2              10.666667         0.3125   \n",
       "107_V5_NS_A1             5.2              10.666667         0.3125   \n",
       "108_V4_NS_A1             0.5               1.800000         0.3125   \n",
       "\n",
       "              protectNorm_PRN protectNorm_FHA geommean_protectNorm  \\\n",
       "SampleID                                                             \n",
       "106_V5_NS_A1           0.3125          1.3750             1.140388   \n",
       "107_V2_NS_A1           1.1250          0.3750             1.783418   \n",
       "107_V3_NS_A1           1.1250          0.3750             1.783418   \n",
       "107_V5_NS_A1           1.1250          0.3750             1.783418   \n",
       "108_V4_NS_A1           0.3125          0.1875             0.449420   \n",
       "\n",
       "              VR_group_v2  \n",
       "SampleID                   \n",
       "106_V5_NS_A1          NVR  \n",
       "107_V2_NS_A1          NVR  \n",
       "107_V3_NS_A1          NVR  \n",
       "107_V5_NS_A1          NVR  \n",
       "108_V4_NS_A1          LVR  \n",
       "\n",
       "[5 rows x 84 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meta = pd.read_csv('../../data/metadata/nasal/nasal_metadata.csv', index_col='SampleID')\n",
    "meta['age_at_collection'] = (pd.to_datetime(meta['CollectionDate']) - pd.to_datetime(meta['DOB'])).dt.days\n",
    "meta = pd.concat([meta,\n",
    "                  pd.read_csv('../../data/metadata/nasal/nasal_abx_usage.csv', index_col='SampleID'),\n",
    "                  pd.read_csv('../../data/metadata/nasal/nasal_titers_yr1.csv', index_col='SampleID')],\n",
    "                 axis=1)\n",
    "meta = meta.loc[~pd.isna(meta['median_mmNorm'])]\n",
    "meta.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "07db70cc-f930-477f-ab3f-27777f8af4f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>101_S1_NS_A1</th>\n",
       "      <th>101_V3_NS_A1</th>\n",
       "      <th>101_V5_NS_A1</th>\n",
       "      <th>102_V1_NS_A1</th>\n",
       "      <th>102_V3_NS_A1</th>\n",
       "      <th>102_V5_NS_A1</th>\n",
       "      <th>102_V6_NS_A1</th>\n",
       "      <th>103_S1_NS_A1</th>\n",
       "      <th>103_S3_NS_A1</th>\n",
       "      <th>103_V10_NS_A1</th>\n",
       "      <th>...</th>\n",
       "      <th>MSA2002_5A</th>\n",
       "      <th>MSA2002_5B</th>\n",
       "      <th>MSA2002_6A</th>\n",
       "      <th>MSA2002_6B</th>\n",
       "      <th>MSA2002_7A</th>\n",
       "      <th>MSA2002_7B</th>\n",
       "      <th>MSA2002_8A</th>\n",
       "      <th>MSA2002_8B</th>\n",
       "      <th>MSA2002_9A</th>\n",
       "      <th>MSA2002_9B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0001</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1593</td>\n",
       "      <td>7320</td>\n",
       "      <td>606</td>\n",
       "      <td>...</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>5845</td>\n",
       "      <td>9876</td>\n",
       "      <td>692</td>\n",
       "      <td>557</td>\n",
       "      <td>783</td>\n",
       "      <td>509</td>\n",
       "      <td>6047</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>114</td>\n",
       "      <td>126</td>\n",
       "      <td>104</td>\n",
       "      <td>115</td>\n",
       "      <td>119</td>\n",
       "      <td>111</td>\n",
       "      <td>168</td>\n",
       "      <td>147</td>\n",
       "      <td>83</td>\n",
       "      <td>103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0003</th>\n",
       "      <td>117</td>\n",
       "      <td>0</td>\n",
       "      <td>879</td>\n",
       "      <td>4392</td>\n",
       "      <td>1428</td>\n",
       "      <td>528</td>\n",
       "      <td>87</td>\n",
       "      <td>877</td>\n",
       "      <td>2642</td>\n",
       "      <td>1498</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0004</th>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>1104</td>\n",
       "      <td>1</td>\n",
       "      <td>6133</td>\n",
       "      <td>475</td>\n",
       "      <td>1</td>\n",
       "      <td>109</td>\n",
       "      <td>2</td>\n",
       "      <td>14</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0005</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4173</td>\n",
       "      <td>24</td>\n",
       "      <td>3121</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 980 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         101_S1_NS_A1  101_V3_NS_A1  101_V5_NS_A1  102_V1_NS_A1  102_V3_NS_A1  \\\n",
       "Otu0001             1             0             0             1             2   \n",
       "Otu0002          5845          9876           692           557           783   \n",
       "Otu0003           117             0           879          4392          1428   \n",
       "Otu0004             9             1          1104             1          6133   \n",
       "Otu0005             0             0             0             0             0   \n",
       "\n",
       "         102_V5_NS_A1  102_V6_NS_A1  103_S1_NS_A1  103_S3_NS_A1  \\\n",
       "Otu0001             0             2          1593          7320   \n",
       "Otu0002           509          6047             1             0   \n",
       "Otu0003           528            87           877          2642   \n",
       "Otu0004           475             1           109             2   \n",
       "Otu0005             1             0          4173            24   \n",
       "\n",
       "         103_V10_NS_A1  ...  MSA2002_5A  MSA2002_5B  MSA2002_6A  MSA2002_6B  \\\n",
       "Otu0001            606  ...           3           4           1           1   \n",
       "Otu0002              3  ...         114         126         104         115   \n",
       "Otu0003           1498  ...           0           0           0           0   \n",
       "Otu0004             14  ...           0           0           1           0   \n",
       "Otu0005           3121  ...           0           0           0           0   \n",
       "\n",
       "         MSA2002_7A  MSA2002_7B  MSA2002_8A  MSA2002_8B  MSA2002_9A  \\\n",
       "Otu0001           0           2           1           1           1   \n",
       "Otu0002         119         111         168         147          83   \n",
       "Otu0003           0           0           0           0           0   \n",
       "Otu0004           0           0           0           0           0   \n",
       "Otu0005           0           0           0           0           0   \n",
       "\n",
       "         MSA2002_9B  \n",
       "Otu0001           0  \n",
       "Otu0002         103  \n",
       "Otu0003           0  \n",
       "Otu0004           0  \n",
       "Otu0005           0  \n",
       "\n",
       "[5 rows x 980 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts = pd.read_csv('../../data/nasal/otu_table.gt10_rar10K.tsv', sep='\\t', index_col=0).transpose()\n",
    "counts.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be1270b2-e08a-4060-a5e3-3008288626f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(775, 84)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/6t/1w2t3qmd1rx81mfw9sq_tfpr0000gn/T/ipykernel_12456/52385233.py:2: FutureWarning: Passing a set as an indexer is deprecated and will raise in a future version. Use a list instead.\n",
      "  meta = meta.loc[in_both].sort_values(['BabyN', 'age_at_collection'])\n"
     ]
    }
   ],
   "source": [
    "in_both = set(meta.index) & set(counts.columns)\n",
    "meta = meta.loc[in_both].sort_values(['BabyN', 'age_at_collection'])\n",
    "print(meta.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "19ccfe3f-9562-46f7-8adf-91892eadfe82",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_v5 = meta.query(\"VisitCode == 'V5'\")\n",
    "counts_v5 = counts[meta_v5.index]\n",
    "counts_v5 = counts_v5.loc[(counts_v5 > 0).sum(axis=1) > counts_v5.shape[1]*.2]\n",
    "\n",
    "meta_v6 = meta.query(\"VisitCode == 'V6'\")\n",
    "counts_v6 = counts[meta_v6.index]\n",
    "counts_v6 = counts_v6.loc[(counts_v6 > 0).sum(axis=1) > counts_v6.shape[1]*.2]\n",
    "\n",
    "meta_v7 = meta.query(\"VisitCode == 'V7'\")\n",
    "counts_v7 = counts[meta_v7.index]\n",
    "counts_v7 = counts_v7.loc[(counts_v7 > 0).sum(axis=1) > counts_v7.shape[1]*.2]\n",
    "\n",
    "meta_v9 = meta.query(\"VisitCode == 'V9'\")\n",
    "counts_v9 = counts[meta_v9.index]\n",
    "counts_v9 = counts_v9.loc[(counts_v9 > 0).sum(axis=1) > counts_v9.shape[1]*.2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7c1b3fe-5a38-4624-b29d-d1fb437276c8",
   "metadata": {},
   "source": [
    "There are some samples which have DTAPHib titers but do not have PCV titers. So here we refilter everything to remove samples that do not have PCV titers for analysis of PCV titers specifically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5c6245d0-c5fa-4693-bab7-d570d30938ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_PCV = meta.loc[~pd.isna(meta['median_mmNorm_PCV'])]\n",
    "meta_PCV.head()\n",
    "\n",
    "meta_PCV_v5 = meta_PCV.query(\"VisitCode == 'V5'\")\n",
    "counts_PCV_v5 = counts[meta_PCV_v5.index]\n",
    "counts_PCV_v5 = counts_PCV_v5.loc[(counts_PCV_v5 > 0).sum(axis=1) > counts_PCV_v5.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v6 = meta_PCV.query(\"VisitCode == 'V6'\")\n",
    "counts_PCV_v6 = counts[meta_PCV_v6.index]\n",
    "counts_PCV_v6 = counts_PCV_v6.loc[(counts_PCV_v6 > 0).sum(axis=1) > counts_PCV_v6.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v7 = meta_PCV.query(\"VisitCode == 'V7'\")\n",
    "counts_PCV_v7 = counts[meta_PCV_v7.index]\n",
    "counts_PCV_v7 = counts_PCV_v7.loc[(counts_PCV_v7 > 0).sum(axis=1) > counts_PCV_v7.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v9 = meta_PCV.query(\"VisitCode == 'V9'\")\n",
    "counts_PCV_v9 = counts[meta_PCV_v9.index]\n",
    "counts_PCV_v9 = counts_PCV_v9.loc[(counts_PCV_v9 > 0).sum(axis=1) > counts_PCV_v9.shape[1]*.2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b6ca75-7f7f-4ffb-a575-5d399ab746e5",
   "metadata": {},
   "source": [
    "## Correlations with median titer values\n",
    "\n",
    "We will use Spearman's R as our correlation metric and use OTU abundances from the 2 month (V5), 4 month (V6), 6 month (V7) and 1 year (V9) time points. 2, 4 and 6 months are when vaccinations are given and 1 year is when titers were measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "07c9fcf2-97a8-47a9-833d-e3eee803110d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0104</th>\n",
       "      <td>-0.272545</td>\n",
       "      <td>0.028060</td>\n",
       "      <td>0.635003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0023</th>\n",
       "      <td>0.257361</td>\n",
       "      <td>0.038485</td>\n",
       "      <td>0.635003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0011</th>\n",
       "      <td>-0.211253</td>\n",
       "      <td>0.091170</td>\n",
       "      <td>0.675647</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0035</th>\n",
       "      <td>-0.207605</td>\n",
       "      <td>0.097031</td>\n",
       "      <td>0.675647</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0068</th>\n",
       "      <td>-0.204427</td>\n",
       "      <td>0.102371</td>\n",
       "      <td>0.675647</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0104 -0.272545  0.028060  0.635003\n",
       "Otu0023  0.257361  0.038485  0.635003\n",
       "Otu0011 -0.211253  0.091170  0.675647\n",
       "Otu0035 -0.207605  0.097031  0.675647\n",
       "Otu0068 -0.204427  0.102371  0.675647"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_correlations = counts_v5.transpose().apply(spearmanr, b=meta_v5['median_mmNorm']).transpose()\n",
    "v5_correlations.columns = ['rho', 'p_value']\n",
    "v5_correlations['p_adj'] = p_adjust(v5_correlations['p_value'])\n",
    "v5_correlations = v5_correlations.sort_values('p_value')\n",
    "v5_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5d4719b9-ee9c-4fae-9806-cf6641bd571c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.299024</td>\n",
       "      <td>0.017282</td>\n",
       "      <td>0.343967</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0007</th>\n",
       "      <td>-0.297600</td>\n",
       "      <td>0.017848</td>\n",
       "      <td>0.343967</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0140</th>\n",
       "      <td>0.270193</td>\n",
       "      <td>0.032222</td>\n",
       "      <td>0.343967</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0014</th>\n",
       "      <td>0.268400</td>\n",
       "      <td>0.033429</td>\n",
       "      <td>0.343967</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0025</th>\n",
       "      <td>-0.261925</td>\n",
       "      <td>0.038108</td>\n",
       "      <td>0.343967</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0029  0.299024  0.017282  0.343967\n",
       "Otu0007 -0.297600  0.017848  0.343967\n",
       "Otu0140  0.270193  0.032222  0.343967\n",
       "Otu0014  0.268400  0.033429  0.343967\n",
       "Otu0025 -0.261925  0.038108  0.343967"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_correlations = counts_v6.transpose().apply(spearmanr, b=meta_v6['median_mmNorm']).transpose()\n",
    "v6_correlations.columns = ['rho', 'p_value']\n",
    "v6_correlations['p_adj'] = p_adjust(v6_correlations['p_value'])\n",
    "v6_correlations = v6_correlations.sort_values('p_value')\n",
    "v6_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dab2c6b9-30ad-4be3-b858-2114a8308d00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0140</th>\n",
       "      <td>-0.341066</td>\n",
       "      <td>0.005431</td>\n",
       "      <td>0.249816</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0027</th>\n",
       "      <td>-0.235725</td>\n",
       "      <td>0.058713</td>\n",
       "      <td>0.683525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0047</th>\n",
       "      <td>0.215799</td>\n",
       "      <td>0.084261</td>\n",
       "      <td>0.683525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0075</th>\n",
       "      <td>-0.214319</td>\n",
       "      <td>0.086465</td>\n",
       "      <td>0.683525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0049</th>\n",
       "      <td>0.201552</td>\n",
       "      <td>0.107397</td>\n",
       "      <td>0.683525</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0140 -0.341066  0.005431  0.249816\n",
       "Otu0027 -0.235725  0.058713  0.683525\n",
       "Otu0047  0.215799  0.084261  0.683525\n",
       "Otu0075 -0.214319  0.086465  0.683525\n",
       "Otu0049  0.201552  0.107397  0.683525"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v7_correlations = counts_v7.transpose().apply(spearmanr, b=meta_v7['median_mmNorm']).transpose()\n",
    "v7_correlations.columns = ['rho', 'p_value']\n",
    "v7_correlations['p_adj'] = p_adjust(v7_correlations['p_value'])\n",
    "v7_correlations = v7_correlations.sort_values('p_value')\n",
    "v7_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "42e16d24-3221-4620-ad29-da3021ed7d0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0065</th>\n",
       "      <td>0.305348</td>\n",
       "      <td>0.011339</td>\n",
       "      <td>0.424769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0023</th>\n",
       "      <td>0.289933</td>\n",
       "      <td>0.016472</td>\n",
       "      <td>0.424769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0016</th>\n",
       "      <td>0.272606</td>\n",
       "      <td>0.024508</td>\n",
       "      <td>0.424769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0010</th>\n",
       "      <td>0.262082</td>\n",
       "      <td>0.030850</td>\n",
       "      <td>0.424769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0115</th>\n",
       "      <td>0.261411</td>\n",
       "      <td>0.031297</td>\n",
       "      <td>0.424769</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0065  0.305348  0.011339  0.424769\n",
       "Otu0023  0.289933  0.016472  0.424769\n",
       "Otu0016  0.272606  0.024508  0.424769\n",
       "Otu0010  0.262082  0.030850  0.424769\n",
       "Otu0115  0.261411  0.031297  0.424769"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_correlations = counts_v9.transpose().apply(spearmanr, b=meta_v9['median_mmNorm']).transpose()\n",
    "v9_correlations.columns = ['rho', 'p_value']\n",
    "v9_correlations['p_adj'] = p_adjust(v9_correlations['p_value'])\n",
    "v9_correlations = v9_correlations.sort_values('p_value')\n",
    "v9_correlations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "967bc78d-674e-47de-a92f-902130f50e84",
   "metadata": {},
   "source": [
    "Across all test nothing is significant after multiple test correction. OTU 140 is significant raw at V6 and V7 but not strong enough to dig into further. The most raw significant results are at 1 year. Another potential indicator that titer is effected by factors present when it is measured?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895a3234-112b-4ec5-8d2e-52bbb6d57f31",
   "metadata": {},
   "source": [
    "## Correlations with median titer group values\n",
    "\n",
    "Now we will split the titers into DTAPHib and PCV and test for significant correlations separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dfaa593b-4625-4c9a-935e-501ea409e857",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0036</th>\n",
       "      <td>0.269762</td>\n",
       "      <td>0.029769</td>\n",
       "      <td>0.584765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.260939</td>\n",
       "      <td>0.035777</td>\n",
       "      <td>0.584765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0023</th>\n",
       "      <td>0.208451</td>\n",
       "      <td>0.095646</td>\n",
       "      <td>0.584765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0022</th>\n",
       "      <td>0.207865</td>\n",
       "      <td>0.096603</td>\n",
       "      <td>0.584765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0068</th>\n",
       "      <td>-0.205359</td>\n",
       "      <td>0.100782</td>\n",
       "      <td>0.584765</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0036  0.269762  0.029769  0.584765\n",
       "Otu0029  0.260939  0.035777  0.584765\n",
       "Otu0023  0.208451  0.095646  0.584765\n",
       "Otu0022  0.207865  0.096603  0.584765\n",
       "Otu0068 -0.205359  0.100782  0.584765"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_DTAPHib_correlations = counts_v5.transpose().apply(spearmanr, b=meta_v5['median_mmNorm_DTAPHib']).transpose()\n",
    "v5_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v5_DTAPHib_correlations['p_adj'] = p_adjust(v5_DTAPHib_correlations['p_value'])\n",
    "v5_DTAPHib_correlations = v5_DTAPHib_correlations.sort_values('p_value')\n",
    "v5_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2a6ff897-fec6-4a82-a4fc-e2d4de36c2e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0023</th>\n",
       "      <td>0.256910</td>\n",
       "      <td>0.047530</td>\n",
       "      <td>0.629053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0104</th>\n",
       "      <td>-0.243666</td>\n",
       "      <td>0.060640</td>\n",
       "      <td>0.629053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0011</th>\n",
       "      <td>-0.233991</td>\n",
       "      <td>0.071951</td>\n",
       "      <td>0.629053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0035</th>\n",
       "      <td>-0.214262</td>\n",
       "      <td>0.100196</td>\n",
       "      <td>0.629053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0079</th>\n",
       "      <td>-0.201801</td>\n",
       "      <td>0.122053</td>\n",
       "      <td>0.629053</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0023  0.256910  0.047530  0.629053\n",
       "Otu0104 -0.243666  0.060640  0.629053\n",
       "Otu0011 -0.233991  0.071951  0.629053\n",
       "Otu0035 -0.214262  0.100196  0.629053\n",
       "Otu0079 -0.201801  0.122053  0.629053"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_PCV_correlations = counts_PCV_v5.transpose().apply(spearmanr, b=meta_PCV_v5['median_mmNorm_PCV']).transpose()\n",
    "v5_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v5_PCV_correlations['p_adj'] = p_adjust(v5_PCV_correlations['p_value'])\n",
    "v5_PCV_correlations = v5_PCV_correlations.sort_values('p_value')\n",
    "v5_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9f5e043a-fafc-402a-a8bc-de132b2a0821",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0059</th>\n",
       "      <td>-0.332116</td>\n",
       "      <td>0.007832</td>\n",
       "      <td>0.302081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0134</th>\n",
       "      <td>0.317639</td>\n",
       "      <td>0.011188</td>\n",
       "      <td>0.302081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.291052</td>\n",
       "      <td>0.020655</td>\n",
       "      <td>0.371783</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0040</th>\n",
       "      <td>0.247697</td>\n",
       "      <td>0.050314</td>\n",
       "      <td>0.679236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0036</th>\n",
       "      <td>0.234684</td>\n",
       "      <td>0.064113</td>\n",
       "      <td>0.692419</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0059 -0.332116  0.007832  0.302081\n",
       "Otu0134  0.317639  0.011188  0.302081\n",
       "Otu0029  0.291052  0.020655  0.371783\n",
       "Otu0040  0.247697  0.050314  0.679236\n",
       "Otu0036  0.234684  0.064113  0.692419"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_DTAPHib_correlations = counts_v6.transpose().apply(spearmanr, b=meta_v6['median_mmNorm_DTAPHib']).transpose()\n",
    "v6_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v6_DTAPHib_correlations['p_adj'] = p_adjust(v6_DTAPHib_correlations['p_value'])\n",
    "v6_DTAPHib_correlations = v6_DTAPHib_correlations.sort_values('p_value')\n",
    "v6_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "72ccf352-6156-4ce1-9831-76cd256e3d9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.329917</td>\n",
       "      <td>0.011435</td>\n",
       "      <td>0.383007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0007</th>\n",
       "      <td>-0.303263</td>\n",
       "      <td>0.020662</td>\n",
       "      <td>0.383007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0034</th>\n",
       "      <td>0.300091</td>\n",
       "      <td>0.022097</td>\n",
       "      <td>0.383007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0025</th>\n",
       "      <td>-0.249917</td>\n",
       "      <td>0.058488</td>\n",
       "      <td>0.687970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0040</th>\n",
       "      <td>0.237373</td>\n",
       "      <td>0.072787</td>\n",
       "      <td>0.687970</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0029  0.329917  0.011435  0.383007\n",
       "Otu0007 -0.303263  0.020662  0.383007\n",
       "Otu0034  0.300091  0.022097  0.383007\n",
       "Otu0025 -0.249917  0.058488  0.687970\n",
       "Otu0040  0.237373  0.072787  0.687970"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_PCV_correlations = counts_PCV_v6.transpose().apply(spearmanr, b=meta_PCV_v6['median_mmNorm_PCV']).transpose()\n",
    "v6_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v6_PCV_correlations['p_adj'] = p_adjust(v6_PCV_correlations['p_value'])\n",
    "v6_PCV_correlations = v6_PCV_correlations.sort_values('p_value')\n",
    "v6_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9af93a96-e47c-4fc9-86d1-abb5e2db5dd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0005</th>\n",
       "      <td>-0.314334</td>\n",
       "      <td>0.010770</td>\n",
       "      <td>0.495437</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0049</th>\n",
       "      <td>0.215649</td>\n",
       "      <td>0.084484</td>\n",
       "      <td>0.841144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0140</th>\n",
       "      <td>-0.203925</td>\n",
       "      <td>0.103236</td>\n",
       "      <td>0.841144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0015</th>\n",
       "      <td>-0.203734</td>\n",
       "      <td>0.103565</td>\n",
       "      <td>0.841144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0037</th>\n",
       "      <td>-0.203222</td>\n",
       "      <td>0.104454</td>\n",
       "      <td>0.841144</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0005 -0.314334  0.010770  0.495437\n",
       "Otu0049  0.215649  0.084484  0.841144\n",
       "Otu0140 -0.203925  0.103236  0.841144\n",
       "Otu0015 -0.203734  0.103565  0.841144\n",
       "Otu0037 -0.203222  0.104454  0.841144"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v7_DTAPHib_correlations = counts_v7.transpose().apply(spearmanr, b=meta_v7['median_mmNorm_DTAPHib']).transpose()\n",
    "v7_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v7_DTAPHib_correlations['p_adj'] = p_adjust(v7_DTAPHib_correlations['p_value'])\n",
    "v7_DTAPHib_correlations = v7_DTAPHib_correlations.sort_values('p_value')\n",
    "v7_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "db30229f-5604-476f-a92c-960aba714300",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0140</th>\n",
       "      <td>-0.319534</td>\n",
       "      <td>0.012066</td>\n",
       "      <td>0.485300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0010</th>\n",
       "      <td>0.294781</td>\n",
       "      <td>0.021100</td>\n",
       "      <td>0.485300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0049</th>\n",
       "      <td>0.253070</td>\n",
       "      <td>0.049090</td>\n",
       "      <td>0.739043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0047</th>\n",
       "      <td>0.235717</td>\n",
       "      <td>0.067432</td>\n",
       "      <td>0.739043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0075</th>\n",
       "      <td>-0.211476</td>\n",
       "      <td>0.101825</td>\n",
       "      <td>0.739043</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0140 -0.319534  0.012066  0.485300\n",
       "Otu0010  0.294781  0.021100  0.485300\n",
       "Otu0049  0.253070  0.049090  0.739043\n",
       "Otu0047  0.235717  0.067432  0.739043\n",
       "Otu0075 -0.211476  0.101825  0.739043"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "v7_PCV_correlations = counts_PCV_v7.transpose().apply(spearmanr, b=meta_PCV_v7['median_mmNorm_PCV']).transpose()\n",
    "v7_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v7_PCV_correlations['p_adj'] = p_adjust(v7_PCV_correlations['p_value'])\n",
    "v7_PCV_correlations = v7_PCV_correlations.sort_values('p_value')\n",
    "v7_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1ec6b34b-73cf-4121-90ae-3bdd5b157f56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0008</th>\n",
       "      <td>-0.374099</td>\n",
       "      <td>0.001674</td>\n",
       "      <td>0.113845</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0123</th>\n",
       "      <td>0.243936</td>\n",
       "      <td>0.045001</td>\n",
       "      <td>0.840422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0095</th>\n",
       "      <td>-0.235700</td>\n",
       "      <td>0.052996</td>\n",
       "      <td>0.840422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0122</th>\n",
       "      <td>0.235481</td>\n",
       "      <td>0.053223</td>\n",
       "      <td>0.840422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0050</th>\n",
       "      <td>-0.227529</td>\n",
       "      <td>0.062038</td>\n",
       "      <td>0.840422</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0008 -0.374099  0.001674  0.113845\n",
       "Otu0123  0.243936  0.045001  0.840422\n",
       "Otu0095 -0.235700  0.052996  0.840422\n",
       "Otu0122  0.235481  0.053223  0.840422\n",
       "Otu0050 -0.227529  0.062038  0.840422"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_DTAPHib_correlations = counts_v9.transpose().apply(spearmanr, b=meta_v9['median_mmNorm_DTAPHib']).transpose()\n",
    "v9_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v9_DTAPHib_correlations['p_adj'] = p_adjust(v9_DTAPHib_correlations['p_value'])\n",
    "v9_DTAPHib_correlations = v9_DTAPHib_correlations.sort_values('p_value')\n",
    "v9_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7f2cd382-5979-429c-af3a-c6eb08d26b56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0065</th>\n",
       "      <td>0.370817</td>\n",
       "      <td>0.002774</td>\n",
       "      <td>0.191372</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0023</th>\n",
       "      <td>0.323433</td>\n",
       "      <td>0.009719</td>\n",
       "      <td>0.335312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0016</th>\n",
       "      <td>0.286060</td>\n",
       "      <td>0.023039</td>\n",
       "      <td>0.375104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0092</th>\n",
       "      <td>0.276624</td>\n",
       "      <td>0.028186</td>\n",
       "      <td>0.375104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.275770</td>\n",
       "      <td>0.028696</td>\n",
       "      <td>0.375104</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0065  0.370817  0.002774  0.191372\n",
       "Otu0023  0.323433  0.009719  0.335312\n",
       "Otu0016  0.286060  0.023039  0.375104\n",
       "Otu0092  0.276624  0.028186  0.375104\n",
       "Otu0029  0.275770  0.028696  0.375104"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_PCV_correlations = counts_PCV_v9.transpose().apply(spearmanr, b=meta_PCV_v9['median_mmNorm_PCV']).transpose()\n",
    "v9_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v9_PCV_correlations['p_adj'] = p_adjust(v9_PCV_correlations['p_value'])\n",
    "v9_PCV_correlations = v9_PCV_correlations.sort_values('p_value')\n",
    "v9_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f36b9cd-bad9-4983-941e-3b13a77e5df3",
   "metadata": {},
   "source": [
    "Still nothing significant after multiple testing correction. Still most raw significance at 1 year (V9) but nothing consistent enough to get excited about despite the low p-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "512c00ba-d5ff-421b-a9f8-2779994ebc29",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
